#!/usr/bin/perl -w
use strict;
use warnings;

if (@ARGV != 1) {
  print "USAGE: $0 output.fasta\n";
  exit(1);
}

my ($fastafile) = @ARGV;
my @seqs = getFastaFrom($fastafile);

if (@seqs == 0) {
  print "Need at least one sequence in $fastafile.\n";
  exit(1);
}

my $sumDistances = 0;

foreach my $seq (@seqs) {
  # get the coding sequence
  $seq =~ s/^.*?([ACGT]{3}.*?[ACGT]{3}).*$/$1/g;
  # homogenize capitalization
  $seq = uc $seq;
  $sumDistances += avgDistance($seq);
}

my $meanOverallDistance = $sumDistances/@seqs;

printf "Mean Distance: %.2f codon", $meanOverallDistance;
($meanOverallDistance == 1) ? print "\n" : print "s\n";

# given a FASTA file, return an array of all sequences from that file (not 
# including sequence labels)
sub getFastaFrom {
  my ($fastafile) = @_;
  my $qty = 0;
  my @seqs;
  open (INPUT, $fastafile) or die "Can't open $fastafile: $!";

  while( my $line = <INPUT> ) {
    if( $line =~ /^\s*$/ ) {
      next;
    } elsif( $line =~ /^>/ ) {
      # increment the quantity of sequences if a label is identified
      $qty++;
    } else {
      # concatenate the line to the current sequence if it is not a label
      $line =~ s/[^ACGT]//gi;
      $seqs[$qty-1] .= $line;
    }
  }
  return @seqs;
}

# Case 1: codon appears more than once. Score += distances, for each distance

# Case 2: codon appears once. Score += sequence length (rationale: appearing
# once is only just slightly better than appearing twice, once at the 
# beginning and once at the end)

# Score /= number of distances
sub avgDistance { 
  my ($seq) = @_;
  # split the sequence into triplets of characters (codons)
  my @codons = ( $seq =~ m/.{3}/g );
  my $sum = 0; # total distance
  my $count = 0; # quantity of distances found

  # search for codons that are the same, record distances in $sum
  for (my $i=0; $i<@codons; $i++) {
    if ( $codons[$i] eq 'visited' ) {
      # this codon was already accounted for
      next;    
    } elsif ( $i == $#codons ) {
      # last codon is unique
      $sum += @codons;
      $count++;
    }
    
    my $codon=$codons[$i]; 
    my $prev=$i; # last position where we saw the codon
    for (my $j=$i+1; $j<@codons; $j++) {
      if ( $codons[$j] eq 'visited' ) {
        # this codon was already accounted for
        next;
      }

      if ( $codons[$j] eq $codon ) {
        # found two equivalent codons in different positions
        $sum += $j-$prev;
        $count++;
        $prev=$j; # new position to measure distance relative to
        $codons[$j] = 'visited'; # mark the second codon as accounted for
      } elsif ( $j == $#codons && $prev == $i ) {
        # a codon other than the last one is unique
        $sum += @codons;
        $count++;
      }
    }   
  }
  # avoid potential divide-by-zero error
  return ($count == 0) ? 0 : $sum/$count;
}

